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ABSTRACT 

An algebraic geometry method is used to calculate the moments of the electron density of 
states as a function of the energy for lattices in the tight binding approximation. Interpreting 
the moments as the Mellin transform of the density allows writing down a formula for the 
density as an inverse Mellin transform. The method is illustrated by working out the density 
function for the two-dimensional square and honeycomb lattices. 
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The tight binding model is a widely used scheme for studying electronic band structure of solids 
[1]. The model is defined by a Hamiltonian quadratic in the electron creation and destruction oper- 
ators indexed by a set of points in the D-dimensional Euclidean space R D , called sites. The sites 
form a lattice A, taken to model a crystal. The physical picture underpinning the model supposes 
that the electrons are tightly bound to a site but may hop from a given site to its neighbouring ones, 
which, for the purpose of the present discussion, are restricted to the nearest neighbours only, with 
proximity defined with respect to distances measured along lattice paths. Thus each physical sys- 
tem is defined by its specific lattice description. The translation symmetry of the lattice permits 
restricting the quasi-momenta k, that is the variables on the reciprocal lattice A, dual to A, to a 
closed subset of the dual H D . The convex hull of this closed subset is called the Brillouin zone. 
The eigenvalues of the tight binding Hamiltonian are invariant functions defined on the Brillouin 
zone. 

We consider a related variant of the tight binding approximation wherein the energy eigen- 
values of electrons are those of a discrete Laplacian associated with the lattice [2]. The discrete 
Laplacian is defined on complex- valued functions / on R D as 

V/(«)= c a c b f(v + a-b), 

a,b e^cA 

where the set A generates the lattice A. The parametres c are taken to be unity on every site. The 
eigenvalues of the Laplacian for this variant are the square of the energy eigenvalues obtained from 
the usual tight binding model. 

Given a lattice A in H D and the energy eigenvalues of the single electron states, its associated 
Green's function, often referred to as the lattice Green's function, can be evaluated and has found 
diverse applications [3-6]. The electronic density of states (DOS) as well as a host of other physical 
quantities of the crystalline solid can be obtained from the Green's function. For example, the 
DOS of the system can be determined from the imaginary part of the Green's function G as p(e) = 
—(1/ir) lim^_ >0 + I m G(e + irj), e denoting the energy eigenvalue. This follows from the definition 
of the density of states as a sum over delta functions 5(E — E n ) over energies, where E n is an 
electron energy eigenvalue. Various techniques have been developed to determine the density of 
states as it contains important physical information, such as electron conductivity in solids [7-9]. 

We use standard methods of algebraic geometry to determine the moments of the density of 
states for the tight binding model in two dimensions. There is a certain naturality in this formu- 
lation. First, the Laplacian is a natural operator and second, periodic functions of two variables, 
like the single electron eigenvalues obtained here, is a means to define a well-studied object in al- 
gebraic geometry, namely a complex algebraic surface, also called an elliptic curve. Thus all two- 
dimensional lattice systems with energy eigenvalues periodic in both directions represent elliptic 
curves. Let us mention that although we restrict to two-dimensional models only, the technique 
used here generalises to higher dimensions. 

The density of states in this approach can be written solely in terms of the combinatorial data 
of the lattice, without requiring the knowledge of electronic wave functions and sum over delta 
functions. There exist algebraic geometry methods for studying elliptic curves using differential 
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equations, known as Picard-Fuchs equations. The solutions to these equations provide an alterna- 
tive way to describe an elliptic curve. From this algebraic geometry insight the electronic density 
of states can also be related to these solutions of the Picard-Fuchs equation which, in our case, is 
a single second order differential equation [10]. Circuits around the singular points of the Picard- 
Fuchs equation are related to the (co)homological properties of the curve. Indeed, the derivation 
of the Picard-Fuchs equation follows from these topological properties. A surface, such as the 
Brillouin zone, which is doubly periodic, is topologically a torus, with two linearly independent 
closed one-forms that are not exact. Let us recall that a closed one-form on a space is one that 
vanishes when operated on by the differential operator d. It can be written locally, though not 
necessarily globally, as df where / is a function on the space. An exact one-form is one that can 
be globally written as df. Such a form becomes identically zero when operated on by the operator 
d. The dimension of the first cohomology group of the space is the number of linearly independent 
closed but not exact one-forms. It is a topological invariant. In the two-dimensional examples that 
we discuss, there are two linearly independent closed but not exact one-forms present. Thus, if we 
start with an arbitrary local expression for a family of one-forms on the surface and differentiate 
with respect to the family parametre z, then every differentiation produces a new one-form. Thus, 
the first and second derivatives along with the original one make three one-forms. If all of these 
one-forms are further constructed to be closed then we know that there must be a linear relationship 
between these three since the first cohomology group has dimension two. This linear relationship 
is the Picard-Fuchs equation. The procedure of constructing the one-forms on algebraic surfaces 
through differentiation with respect to the family parametre and discarding exact one-forms at each 
step has been used earlier in various contexts [10-12]. 

We illustrate this approach in two examples, namely, the two-dimensional tight binding model 
for the square and honeycomb lattices. The honeycomb case represents graphene which is a system 
of considerable current interest. In these two cases we show the two steps used to determine 
the electron energy density of states. First, algebraic geometry is used to determine the energy 
moments of the density of states which are interpreted as Mellin transforms. Next, we use the 
powerful techniques of inverting Mellin transforms to determine an analytic expression for the 
density of states. Indeed, the advantage of the present approach lies in obtaining the density of 
states as an inverse Mellin transform, simplifying numerical evaluations for any value of the energy. 

Let us start by briefly discussing the general combinatorial set up to fix notation. We shall also 
identify the physical quantities, in particular, the density of states, in terms of the combinatorial 
data. We restrict the discussion to two-dimensional cases. Generalisation to higher dimensions 
may be considered following known results [2]. The model we consider is described by a finite 
subset A of Z 2 . The lattice A is then obtained by taking the Z-span of the difference of points in 
A, that is 

A = Z-span{a - b\ a, b e A}. (1) 

In other words, the set A is obtained as marking one of the lattice points of the model as the 
origin and collecting the points connected to it by a single path in the lattice. For example, the set 
A = {(-1,0), (1, 0), (0, -1), (0, 1)} for the square lattice, while A= {(1, 0), (0, 1), (-1, -1)} 
for the honeycomb lattice in two dimensions. The lattice constant is taken to be unity throughout. 
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On the points of A we consider the distribution given by a sum of Dirac deltas as 

V = <*«■ (2) 
a eA 

Being supported solely on the lattice points, this embodies a crystal in the tight binding approx- 
imation. The delta functions may, in principle, have different weights at different points, but we 
shall not consider that here. The Fourier transform of V is given by 

V(k) = e- 27Tik - a , (3) 
a eA 

where the quasi-momenta k = (k 1: k 2 ) are valued in the reciprocal lattice 

A = {k G R 2 |fc • (a- b) G Z,Va,b G A} (4) 

dual to A. The eigenvalues of the discrete Laplacian V based on A are then written in terms of the 
quasi-momenta as the dispersion relation 

E{k) 2 := \V{k)\ 2 = cos27ik- (a-b). (5) 

a,b eA 

The energy E is periodic with period lattice A thus descending to a function on the Brillouin zone 
U A ~ R 2 /A, which has the topology of a torus. Let us introduce complex variables x, y and define 
a Laurent polynomial [2] 

W(x,y)= xa ~ b i ( 6 ) 

a,beA 

associated to the set A, satisfying \V{k)\ 2 = W(e 2mk \ e 2mh ' 2 ), where x = (xi,x 2 ) = (x,y) and 
2 2 , for A G A. The number of states, denoted V(e), is given by the normalized volume 
of the Brillouin zone such that \V(k)\ 2 < e. Let us remark that, as mentioned before, by equation 
(5), the parametre e is the square of the energy obtained from an usual tight binding model. The 
Hilbert transform of the differential dV is defined as the integral of the resolvent l/(z — e) with 
respect to the measure defined by dV over the real line as 



H(z) = I 



z — e 

1 f 1 dx dy 



(2iri) 2 J z — W(x,y) x y 



(8) 



where z is a complex parametre. 

The function H(z) in (8) is the period of a differential one-form along a one-cycle on the 
hypersurface given by z = W(x, y) in (C*) 2 . It is obtained as a solution to a Picard-Fuchs equation 
in the form of a Laurent series in the complex variable z which, according to (7), is given in terms 
of moments a n as 

oo 

H(z) = J2^nZ- 1 - n . (9) 

?1=0 
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The moments can be calculated by either substituting (9) in the Picard-Fuchs equation or by using 
the residue theorem in (8) as 



a n = constant term of the Laurent polynomial W(x, y) n . 

From the moments one can calculate the lattice Green's function and hence the density of states [7, 
8]. Although these methods do not yield explicit formulas, they lead to systematic approximation 
schemes that can be numerically implemented in an efficient manner. We shall consider a different 
way to obtain the density of states from H(z) which yields explicit formulas. The idea is to expand 
(7) in a geometric series in ej z as 

n=0 ' 



R 



oo „ 
n=0 ^ R 



where we defined the density of states p(e) = dV/de. Comparing with (9) we conclude, 

a n = [ p(e)e n de. (12) 
Jr. 

We now make our simple but important observation, namely, the moments a n of the density of 
states p(e) can be interpreted as the Mellin transform of ep(e) if we replace the integers n by a 
complex variable s. An immediate consequence of this remark is, as emphasised before, that an 
expression for the density of states can be easily written down as the inverse Mellin transform of 
a n = a(s). We have the formula 

co+ioo 
co — ioo 

where the line integral is evaluated along a vertical line in the complex plane and c is an appropri- 
ate real constant. This approach thus gives an explicit formula for the density of states in terms of a 
function determined by the methods of algebraic geometry. Moreover, it allows us to calculate the 
density of states for any value of e, large or small, by choosing appropriate contours in the s-plane. 
In order to use this method we need to be able to replace the discrete set a n by a function a(s) of a 
complex variable s. For the cases that we study there is a natural way of doing this. We shall now 
consider two examples. 



Example 1: Square lattice 

For the square lattice the set of generating points in Z 2 is A = {(—1, 0), (1, 0), (0, —1), (0, 1)}. 
This corresponds to the polynomial 

W = (x + l/x + y+l/y) 2 (14) 
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in the coordinate ring C[x,x 1 ,y,y x ] . Then the dispersion relation is obtained to be 

E(k) 2 = 4 + 2 cos 1vk x + 2 cos 2nk 2 . 



(15) 



We shall evaluate the resolvent H(z) defined in (8). Writing the complex variables x, y in terms 
of the homogeneous coordinates of a two-dimensional complex projective space P 2 as x = Vi/vq 

and y — v 2 /vo, vq ^ 1, we rewrite H as 



H(z) 



V ViV 2 Q 



z(v viv 2 ) 2 - (vi + v 2 ) 2 {v 1 v 2 + vl) 2 ' 

_J_(f 5 + f 

2z \J (v viv 2 ) - t(Vi + v 2 )(viv 2 + vl) J 



n 



{v VlV 2 ) -t(vi+ V 2 )(V!V 2 + vl) J (V V!V 2 ) + t(vi + V 2 ) (v X V 2 + Vq) J ' 

(16) 

where f2 = v dv 1 A dv 2 — Vidv A dt> 2 + v 2 dv A (if! is the canonical 2-form on P 2 and we defined 
t — 1/ ^/i. The Picard-Fuchs equations of both the varieties 



t(vi + V 2 )(viV 2 + 1>q) ± V V!V 2 = 0, 



are the same, namely, 



d 2 w 1 - 48t 2 cfcj 16c7 



+ 



+ 



0. 



(17) 



(18) 



dt 2 t - 16t 3 dt lQt 2 - 1 
Thus, series solutions with only the terms with even powers of t survive. This Fuchsian equation 
has two solutions which can be obtained as series in t by the Frobenius' method. The two solutions 
are 



w^t) = 2 Fi(l/2, 1/2; 1; 16£ 2 ), 

d 



w 2 {t) = (logt + 2 1(^2)^(1/2, 1/2; 1; 16t 2 ) + -L ^ _ 

n=0 



r(a + n + 1/2)' 
r(a + m + l) 2 



Q = 



(19) 

;i6t 2 ) n . 

(20) 



Here 2 -Fi( a i> 02; &i! #) represents a hypergeometric function defined by the series 

tp r 7 \ (01)71(02)0 n 
2 F 1 (a 1 , a 2 ; b 1 ;x)=2_^ -77^ v~x , 



n=0 



(21) 



where (a) n = a(a + I). ..(a + n — 1) is the Pochhammer symbol. The two integrals in H are then 
linear combinations of the two solutions, namely, 



H(z) = fVl(l/2, 1/2; 1; 



+ | (log2 _ 41og2)2Fl(1/2 , 1/2;1; ^_^_ E; |_ 

n=0 



r(a + n + 1/2) 2 
T(a + n + l) 2 



a=0 



,16, n 

(22) 



5 



where c\ and c 2 are arbitrary constants. Instead of trying to determine the constants from boundary 
conditions, we shall recourse to the calculation of moments to determine the density of states. This 
entails direct evaluation of the integral (8) using residues. Since the constant term in the expansion 
of W n is ( 2 ") 2 , we have 

{2n\ 2 



n / 

T(2n + l) 2 

r(i + n) 4 (23) 

(2n) 2 T(2n) 2 



n 2 T(n) 2 T(l + n) 2 
i r(n + l/2) 2 
7T T(l + n) 2 ' 

where we used the duplication formula n 1 / 2 Y(2x) = 2 2x ~ l T(x)T(x + 1/2) in the last step. By (9) 
this gives the resolvent H as 

oo 

n=0 (24) 

= 1^(1/2,1/2;!;^). 

Zz z 

By (13), the density of states is then obtained as the inverse Mellin transform 

co+ioo 

p(e) = —- I F(g + V2)2 16'e- 1 - (25) 
P[ ' 2m tt J r(l + s) 2 ' 1 } 

eg— ioo 

where the integrand is derived from (23) by substituting s for n. Choosing c = and closing the 



20 - 




Figure 1: Behaviour of density of states p{e) near e = for the square lattice 



contour with a semicircular arc on the left so as to obtain an expression valid near e = 0, we get 

, , 4 log 2 — log e . . , . . 7r *^ d 

= 4v? 2fl(1/2 - 1/2; 1; f / 16 > + 4^^ z 

v v n=0 

(26) 



r(i + s) 2 r(i/2-s) 2 



.s=-l /2-n 



e n " 
16, 
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whose leading behaviour near e = is shown in Figure 1. Let us point out that there is a lower 
limit to the admissible range of energy e determined by the sample geometry. The density of states 
plotted integrated over the whole range of energy does not depend on this physical cut off since 
the integral is finite even without a cutoff. The same holds good for the honeycomb lattice as well 
to which we now turn as our next example. 

Example 2: honeycomb lattice 

For the honeycomb lattice we have A = {(1,0), (0, 1), (—1, —1)}, leading to the Laurent polyno- 
mial 

W(x,y) = (x + y + —)(- + - + xy) (27) 
xy x y 

in C[x, x~ l ,y, y' 1 ]. The dispersion relation is 

E(k) 2 = 3 + cos(2tt(A;i - fc 2 )) + cos(2tt(2A;i + k 2 )) + cos(27r(/ci + 2/c 2 )) , (28) 
which upon a change of basis of the reciprocal lattice 

fci = (VSk x + 3f%)/6, k 2 — (VSk x - 3/%)/6, (29) 
yields the more usual form [13] 

E(nf = 1 + 4 cos 2 TiKy + 4 cos nK y cos ttV3k x . (30) 
Let us define the homogeneous coordinates of a P 2 , namely [v : v\ : v 2 ], related to x, y by 

x 2 y = v 1 /v , xy 2 = v /v 2 . (31) 
Substituting these in (8) we obtain the resolvent 

H(z)= [ -, (32) 

J ZVqViV 2 - (v + Vi + W 2 )(^0^1 + ViV 2 + V 2 Vo) 

solving the Picard-Fuchs equation 

d 2 zu 9-20^ + 3z 2 dzu (z - 3)zu 



dz 2 z(9 - Wz + z 2 ) dz z(9 - lOz + z 2 ) 



0. (33) 



Again, instead of writing down all the solutions of this equation, it suffices for our purposes to 
consider the moments. The constant term in the expansion of W gives the moments [2] 



n\ 2 ( 2j 



3) \J 



3=0 

r(i + n) 2 r(2j + i) 



. • //-./) 2 i'ii • ./)' 

1 ^ r(i + n) 2 r(j + i/2) 



v^^ r (i + n-j) 2 r(i + J ) 3 



4- 7 
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1 y 1 1 



r(i + n) 2 r(j + i/2) 



v^^ r (i + ^-j) 2 r(i +J ) 3 



—A 1 



(34) 



= 3 F 2 (l/2,-ra,-ra;l,l;4), (35) 
where 3 .F 2 (ai, a 2 , a 3 ; b\, b 2 ; x) is the generalised hypergeometric function defined by the series 



jF 2 (ai, a 2 , a 3 ; £>i, b 2 ; x) = > ( 77 v 



ai)n(a 2 ) n (a 3 ) n n 



'2)n n\ 



X 



(36) 



The duplication formula has been used in deriving the expression (34) and the sum has been ex- 
tended to all integral values of j since 1/T(1 + n — j) vanishes for all j > n + 1. As before, 
equation (33) is solved with 



(37) 



n=0 



Then the density of states is expressed in terms of the inverse Mellin transform of a(s) as 



1 1 



7r 2ni 



Cn+lOO 
OO „ , . , x 2 



k=0„ 



r(l + s) 2 r(fc + l/2) 

r(i + s - k) 2 r(i + ky 



k-l-s 



-4 K e 



(38) 



cq— ICO 



We can also write the sum over k as an integral, as 

co+ioo 



\f Ti 2tii 



dt 



d r(i + 8 )*r(t + i/2)r(-fl 



(39) 



r(i + s -t) 2 r(i + kf 

C co— too 

where the contour C is chosen so as to go parallel to the imaginary axis and closing on the right to 
enclose integers t — 0, 1, 2, • • • on the t -plane. Now reversing the order of the integrations we first 
work evaluate the integral over t by closing the contour on the left, Re(£) < 0, so that we pick up 
contributions from the poles of T(l/2 + t) at t — — 1/2 — k, for positive integers k. This leads to 



00 

Ly 



CQ+lOO 



2v/7r27ri 



fc=0 



r(l + s) 2 r(fc + l/2) 



T(3/2 + fc + s) 2 r(l/2 - fc) 2 r(l + fc) 



4" fc e 



k-l-s 



(40) 



Cq— JOO 



In order to derive a power series in e, we note that T(l + s) 2 has double poles at s = — 1 — n, for 
positive integral n. Thus, performing the integral by closing the contour on the left we obtain 



v ra=0 n=0 



r(fc + l/2)4" fc e 



fc -1-s 



r(3/2 + fc + s ) 2 r(i/2 - fc) 2 r(i + fc)r(-s) 



(41) 



s= — 1— n 



which can be rewritten as 

1 00 00 



EE 



(l/4) fc r(l/2 + fc) e „ 



V^^r(i/2-*)'r(i + *) 



r(l/2 + fc-n) 2 r(l + n) 2 



loge 



+ 



d_ 

ds 



r(3/2 + fc + s ) 2 r(-s) 2 



(42) 



s=— 1— ?v 
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Figure 2: Behaviour of density of states p(e) near e = for the honeycomb lattice 

The behaviour of the density of states near e = is plotted in Figure 2. 

To summarise, we have discussed the density of states of the two-dimensional nearest neigh- 
bour tight binding Hamiltonian from an algebraic geometry viewpoint. We have discussed two 
examples based on the two-dimensional square and honeycomb lattices. The density of states is 
obtained as a function of energy. The Hilbert transform of the density of states is the resolvent that 
satisfies Picard-Fuchs equations of algebraic varieties that correspond to the lattices in a combina- 
torial fashion. Explicit expressions are given for small energies in terms of infinite series, involving 
hypergeometric functions. Let us note that the Brillouin zone corresponding to each of the lattices 
is a topological torus. Thus qualitative features of the results may be understood in topological 
terms. We intend to present details of these topological arguments in a future work. A practical 
advantage of this approach is that it allows evaluation of density of states in any domain of energy 
by appropriate choice of contours suitably in the integrals. The resulting infinite series obtained 
converge rather fast and may be easily evaluated numerically. Finally, let us mention that we have 
presented our calculations in the context of electrons in a crystalline medium but the results ob- 
tained are also applicable for the density of states of a system of phonons where a tight binding 
nearest neighbour model is appropriate [14]. 
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